setwd("C://Users/hpeter.INTRANET/Desktop/multivariate stats R/2024/scripts/")


# CLASSIFICATION AND REGRESSION TREES -----------------------------------------------------


# example classification tree using the iris data set -----------------------------------------------


library(partykit)
?ctree
vignette("ctree", package = "partykit")

data(iris)
head(iris)


irisct <- ctree(Species ~ Petal.Length + Petal.Width + Sepal.Length + Sepal.Width, data=iris)
irisct
plot(irisct)

irisct2 <- ctree(Species ~ .,data = iris)
plot(irisct2)
plot(irisct2, inner_panel = node_barplot)


# confusion matrix
table(predict(irisct), iris$Species)

#predict for new iris specimen:
new.iris1<-data.frame(Sepal.Length=4.7,Sepal.Width=3.4,Petal.Length=1.8,Petal.Width=0.2)
predict(irisct, newdata=new.iris1)

new.iris2<-data.frame(Sepal.Length=5.7,Sepal.Width=2.6,Petal.Length=4.1,Petal.Width=1.3)
predict(irisct, newdata=new.iris2)



# example Doubs fish community --------------------------------------------

library(vegan)
load("Doubs.RData")

# Remove empty site 8
spe <- spe[-8, ]
env <- env[-8, ]
spa <- spa[-8, ]
latlong <- latlong[-8, ]


# Compute matrix of chord distance among sites
spe.norm <- decostand(spe, "normalize")
spe.ch <- vegdist(spe.norm, "euc")

# Cluster sites based on fish species abundance using Ward's method
spe.cw <- hclust(spe.ch, method = "ward.D2")
plot(spe.cw)

# partition dendrogram into groups
?cutree
spe.cw.g <- cutree(spe.cw, 4)
spe.cw.g

library(dendextend)
colors_to_use <- spe.cw.g      #define colors and sort according to tips in dendrogram
colors_to_use<-colors_to_use[order.dendrogram(as.dendrogram(spe.cw))]
spe.cw<-as.dendrogram(spe.cw)
labels_colors(spe.cw) <- colors_to_use   #change color of tip labels depending on best partition
plot(spe.cw)


#combining all data
all.dat<-cbind(spe.cw.g,env,spe, spa)
head(all.dat)

#spatial distribution of sample groups
plot(spa, type="l",col="lightblue",lwd=3)
points(all.dat$X,all.dat$Y,pch=16, cex=2.5, col=all.dat$spe.cw.g, xlab="X coordinate", ylab="Y coordinate")

#use some environmental parameter to predict groups
class.groups=ctree(as.factor(spe.cw.g)~ ele+slo+dis+pH+har+pho+nit+amm+oxy, data = all.dat)
plot(class.groups)

plot(spa, type="l",col="lightblue",lwd=3)
points(all.dat$X,all.dat$Y, pch=16, cex=all.dat$amm*10, col=all.dat$spe.cw.g, xlab="X coordinate", ylab="Y coordinate")


#apply the same logic to generate a regression tree (quantitative response) for BOD 
class.bod<-ctree(bod ~ ele + slo + dis +pH +har+pho +nit + amm +oxy, data=all.dat)
plot(class.bod)

# estimate fish species richness (quantitative) and use a 
# regression tree

all.dat$spec.rich<-rowSums(spe>0)
?diversity #Shannon, Simpson 
vignette("diversity-vegan",package = "vegan")
all.dat$H<-diversity(spe)
head(all.dat)

#TASK
# build regression trees and inspect results for species richness and Shannon H.




# RANDOM FOREST CLASSIFIER ------------------------------------------------
library(randomForest)
?randomForest

env1=env[,-1] # remove distance from source

rf.envcw = randomForest(as.factor(spe.cw.g)~., env1, ntree=500, mtry=3, importance=TRUE,
                        na.action=na.omit,do.trace=100,proximity=T)
rf.envcw


# Variable importance
importance(rf.envcw)
varImpPlot(rf.envcw)


# Partial dependence plots
par(mfrow=c(2,3))
partialPlot(rf.envcw, env1, ele)
partialPlot(rf.envcw, env1, dis)
partialPlot(rf.envcw, env1, amm)
partialPlot(rf.envcw, env1, pho)
partialPlot(rf.envcw, env1, oxy)
partialPlot(rf.envcw, env1, nit)


#TASK 2: 
# build RF using BOD (Biological Oxygen Demand) and the remaining environmental variables
# Plot and interpret the results of the randomForest object (VarImPlot, Partial dependence).